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A New Exact Method for Dynamical Fermion Computations with 
Non-Local Actions 

Ivan Horvath a and A. D. Kennedy b * t and Stefan Sint b * 
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We introduce a new algorithm which we call the Rational Hybrid Monte Carlo Algorithm (RHMC) . This method 
uses a rational approximation to the fermionic kernel together with a noisy (Kennedy-Kuti Jl]) acceptance step 
to give an efficient algorithm with no molecular dynamics integration step-size errors. 
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1. Introduction 

The RHMC algorithm is an exact method for 
generating lattice configurations distributed ac- 
cording to some non-local action. Both the terms 
'non-local' and 'exact' require some explanation. 

By a non-local action we mean one which is 
constructed from some continuous function of a 
local kernel, for example ip^/Mip where M. is 
the usual staggered kernel; we do not mean an 
action with complicated long-range interactions 
and many coupling constants such as a 'perfect' 
action. Two examples of interesting non-local ac- 
tions are two flavours of staggered fermions, and 
solutions of the Ginsparg- Wilson relation such as 
Neuberger's action 

Our goal is not to consider the virtues or other- 
wise of such actions, but only to discuss how they 
might be efficiently simulated. 

'Exact' algorithms in our sense are ones which 
do not require a zero integration step-size extrap- 
olation. This is important in the case of two 
flavours of staggered fermions in order to deter- 
mine the order of phase transitions, for example. 
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2. Ideas underlying the algorithm 

Our method is a Hybrid Monte Carlo (HMC) 
algorithm with two new ingredients, a cheap ac- 
cept/reject step, and a cheap force computation. 

For the accept/reject step, which is needed to 
make the algorithm exact, the problem is that it 
is too expensive to compute the change in energy 
SH exactly. We therefore apply the Kennedy- 
Kuti method [El and use a linear accept/reject 
step by defining an ordering of the fields indepen- 
dent of the estimator of the action. The number 
of violations of < P acc < 1 can be counted ex- 
plicitly and easily be made negligibly small by a 
suitable choice of algorithmic parameters. In or- 
der to produce an unbiased stochastic estimate of 
cxp tr In M. we sum the series expansions for e x 
and ln(l — x) stochastically [f|. 

2.1. Noisy Force 

The computation of the force needed for the 
molecular dynamics (MD) evolution is also pro- 
hibitively expensive for non-local actions, and we 
first investigated modifying the Hybrid Molecu- 
lar Dynamics (HMD) i?o algorithm by adding an 
accept/reject step to make it exact. This did not 
lead to a feasible algorithm because the integra- 
tion errors grow as 5H cx VSt, and thus the cost 
grows as V(VSt)(^/8t) = V 2 £ where St is the 
integration step size, £ is the correlation length 
in units of MD time and V is the lattice volume. 
This cannot be improved using higher-order in- 
tegration schemes because the errors are intrinsi- 
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cally due to the noise. The R algorithm g is not 
area-preserving or reversible so cannot be made 
exact in any obvious fashion. 

2.2. Polynomial Force 

Instead of using a noisy force we then observed 
that the HMC algorithm does not require that 
we integrate the classical equations of motion for 
the Hamiltonian corresponding to the action we 
are interested in: any area-preserving reversible 
mapping on phase space suffices to give a valid 
algorithm. Using an MD Hamiltonian which ap- 
proximates the desired one sufficiently well to give 
a good acceptance rate is a good choice, and fol- 
lowing Luscher j?J we can do this using a polyno- 
mial in M with a minimax error over the compact 
spectrum of M. which falls exponentially with the 
degree of the polynomial. 

The existence of such optimal polynomials was 
shown by Chebyshev, and a truncated expansion 
in Chebyshev polynomials is usually very good, 
although in general not optimal. The optimal 
polynomials can be computed easily using the Re- 
mez algorithm. It is not obvious that the mini- 
max norm is the best for our purposes, indeed 
Montvay uses an L 2 norm instead, but it seems 
to be a safe and adequate choice. 

Such polynomial actions have been used before 
p|-^0[|, but not in combination with a noisy ac- 
ceptance step for non-local actions. The major 
difficulty with using such high-degree polynomial 
actions is that they can be very sensitive to round- 
ing errors in numerical computation, and much 
work has gone into choosing orderings of the roots 
of the polynomial which minimise these difficul- 
ties Q. 

2.3. Rational Force 

It is well-known in the numerical analysis liter- 
ature that rational function approximations not 
only share the exponential convergence property 
of polynomial approximations, but also usually 
reach a sufficiently small minimax error for much 
lower degree. The existence of optimal rational 
approximations is easily established using Cheby- 
shev's arguments, but in this case we must use 
an iterative method such as the Remez algorithm 
to find the optimal approximation. If we choose 
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Figure I . Comparison of minimax errors for ratio- 
nal functions of degree [n,d\. The dashed curves 
(with d — 0) are polynomials. 



to use a diagonal rational approximation, that 
is one in which the numerator and denominator 
are polynomials of equal degree, then it turns out 
that the roots of both are all real and negative. 
An example of such an approximation is 



•fx 



=S 0.3904603901 X 

(z+2. 3475661045) {x+0. 1048344600) (x+0. 0073063814) 
( rr+0. 4105999719) (x+0. 0286165446) (rc+0. 0012779193) ' 



The evaluation of such a rational approximation 
for a matrix requires several conjugate gradient 
inversions with increasing masses.^] Figure [l] il- 
lustrates the quality of the minimax approxima- 
tions. 

3. Numerical results 

The results of our numerical tests of the RHMC 
algorithm are shown in Figure || for the case of 
four flavours of staggered fermions, where the re- 
sults may be compared with a conventional HMC 
computation, and in Figure || for the case of two 
staggered flavours, where there is no other exact 



5 We are also investigating using a partial fraction expan- 
sion of the rational function and using a multiple mass 
solver. 
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Figure 2. Numerical Results for Nf = 4. 



Figure 3. Numerical Results for Nf = 2. 



method available for comparison. 

Each of these figures is composed of four 
graphs: the top one shows the mean acceptance 
rate, which for our choice of parameters should 
be exactly 70%; the second one shows the num- 
ber of violations, when this is zero the method is 
exact; the third one shows the mean plaquette; 
and the bottom one the largest eigenvalue of the 
staggered Dirac operator. The solid circles are 
results using a [3, 3] rational approximation for 
the force and the open ones are for a [32, 0] poly- 
nomial approximation. The lines are the results 
from a conventional HMC computation when this 
was possible. 

Results are shown both for the algorithm with 
the noisy accept/reject step (indicated by the suf- 
fix N) and for an inexact version without the ac- 
cept/reject step (indicated by the suffix D). The 
St 2 step size errors in the latter are evident, as 
are the difficulties in carrying out an accurate zero 
step size extrapolation. 
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